Magnetic susceptibility of a Cu02 plane in the La2Cu04 system: 
I. RPA treatment of the Dzyaloshinskii-Moriya Interactions 

K. V. TabunshchylQ 

Department of Physics, Queen's University, Kingston ON K7L 3N6 Canada and 
Institute for Condensed Matter Physics, Lviv, Ukraine 

R. J. Gooding 

Department of Physics, Queen's University, Kingston ON K7L 3N6 Canada 
(Dated: February 2, 2008) 

Motivated by recent experiments on undoped La2Cu04, which found pronounced temperature- 
dependent anisotropies in the low-field magnetic susceptibility, we have investigated a two- 
dimensional square lattice of 5* = 1/2 spins that interact via Heisenberg exchange plus the symmetric 
and anti-symmetric Dzyaloshinskii-Moriya anisotropies. We describe the transition to a state with 
long-ranged order, and find the spin-wave excitations, with a mean-field theory, linear spin-wave 
analysis, and using Tyablikov's RPA decoupling scheme. We find the different components of the 
susceptibility within all of these approximations, both below and above the Neel temperature, and 
obtain evidence of strong quantum fiuctuations and spin-wave interactions in a broad temperature 
region near the transition. 



I. INTRODUCTION 

Quantum magnetism of low-dimensional systems has 
attracted considerable attention in recent years, in part 
due to the strong interest in the cuprate superconduc- 
tors. For example, it has been postulated that a strong 
antiferromagnetic (AF) exchange interaction may be re- 
sponsible for the high-temperature superconductivity in 
these compounds i 

The ubiquitous structural and electronic constituent of 
this latter class of materials is the Cu02 plane, and in 
this paper we consider the magnetic properties of such 
planes in their undoped state. In particular, we con- 
sider the temperature dependence of the static, uniform 
magnetic susceptibility for a single plane in an undoped 
La2Cu04 crystal. This system is known to be an AF 
insulator with a very simple structure, namely it can be 
approximately thought of as one Cu02 plane stacked be- 
tween LaO planes, with this structural unit repeated, in 
a body-centred tetragonal pattern, throughout all space. 
However, a small orthorhombic distortion introduces im- 
portant spin-orbit couplings into the magnetic Hamilto- 
nian, leading to an AF state with a weak canted fer- 
romagnetic moment. These spin-orbit interactions are 
central to the results presented in this paper. 

As was known from the start of research on the cuprate 
superconductors, a complete knowledge of the proper- 
ties of the spin- 1 quantum Heisenberg AF on a square 
lattice is an absolute necessity^ However, some experi- 
ments have demonstrated that a complete description of 
the magnetic behaviour found in, e.g. La2Cu04, requires 
additional physics. Examples include (i) weak ferro- 
magnctism in the low-temperature orthorhombic (LTO) 
phasc(^ (ii) spin wave gaps with in- and out-of-plane 
modesj^ and perhaps most importantly, (iii) the unusual 
anisotropy of the magnetic susceptibility observed by 
Lavrov, Ando, Komiya and Tsukada»& It was this latter 



experiment that led us to complete a sequence of the- 
oretical investigations on a model that should describe 
such a three-dimensional array of such Cu02 planes mod- 
elling La2Cu04, a structure similar to those found in 
many cuprate superconductors. This manuscript sum- 
marizes the first of these studies, that concerned with 
a single Cu02 plane, with this plane described by a 
near-neighbour Heisenberg model plus spin-orbit cou- 
plings as embodied by the antisymmetric and symmetric 
Dzyaloshinskii-Moriya (DM) interactionsi^ 

An important point needs to be raised to clarify the 
applicability of this work to a real physical system, such 
as La2Cu04. Firstly, note that according to the Mcrmin- 
Wagner theorem a two-dimensional (2D) system with a 
continuous symmetry cannot undergo a continuous phase 
transition, at any nonzero temperature, to a state with 
true long-ranged order. However, when one includes both 
the antisymmetric and symmetric DM interactions this 
symmetry is lifted, and thus the model that we study in 
this paper will have a true phase transition to an ordered 
phase at some nonzero temperature, which we shall la- 
bel by T/v, in analogy to the Neel ordering temperature 
of a pure antiferromagnet. So, the ordered phase for our 
model of a single plane will include a weak ferromagnetic 
canted moment, as well as long-ranged AF order. Note 
that current estimates^ of another interaction present in 
the physical La2Cu04 system, that being a very weak 
AF interlayer coupling which is usually denoted by J_l, 
is that this energy scale is close to that of the DM in- 
teractions, and thus it is likely that both this exchange 
and the DM interactions are roughly equally responsi- 
ble for the observed transition. This serves to emphasize 
that our study of a single plane is not expected to accu- 
rately explain all of the observed magnetic properties of 
La2Cu04; in fact, this work stands alone as a theoretical 
study of an isolated plane, but it is of considerable inter- 
est to learn which experimental data can and which data 
can not be explained by such a single-plane model. 
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We focus on the role of the DM interaction be- 
tween the neighbouring spins in a Cu02 plane. This 
interaction arises from the orthorhombic distortion in 
La2Cu04 (which is associated with the small tilt of the 
CuOg octahedra) together with the spin-orbit interac- 
tion. The DM interaction leads to a small canting of 
the Cu spins out of the plane, so that the weak ferro- 
magnetic order appears in each Cu02 plane, and subse- 
quently allows for the formation of 3D AF order. This 
allows one to observe a pronounced peak in zero-field 
magnetic susceptibility^ x'^(r), and the earliest work on 
the importance of this interaction focused on DM physics. 
To be specific, Thio et almii^ analyzed their susceptibility 
data using a Landau theory expanded to sixth order, for 
a 2D Heisenberg antiferromagnet with interlayer coupling 
and the DM-generated terms. They obtained reasonable 
fits of their theory to the susceptibility and field depen- 
dent magnetization data, and deduced parameters which 
characterized magnetic properties of the La2Cu04 sys- 
tem. As we shall explain below, we believe that the 
necessity of incorporating such higher-order terms into 
their fits is suggestive of the important role played by 
spin-wave interactions, a conclusion consistent with the 
results presented in this, and our future manuscripts on 
this problem. 

Investigations of the magnetic ground state of 
La2Cu04 were performed by several groups of authors, 
usually within the framework of the linear spin-wave 
(SW) theory. The calculations were based on an effec- 
tive model Hamiltonian derived by Moriya's perturbation 
theory^ applied to Hubbard type Hamiltonians by taking 
into account the spin-orbit coupling. In the most gen- 
eral form, the effective spin Hamiltonian, in addition to 
the isotropic exchange interaction, includes the above- 
mentioned antisymmetric and symmetric DM interac- 
tions. The first microscopic derivation of the spin Hamil- 
tonian was performed by Coffey, Rice, and Zhang^i^ they 
estimated the antisymmetric DM coupling constants and 
showed that when the DM vectors alternate a net fer- 
romagnetic moment may be generated in the ground 
state. Shekhtman, Entin-Wohlman, and Aharony^S sub- 
sequently showed that the symmetric anisotropics con- 
tribute to the magnetic energy in the same order as 
the antisymmetric DM anisotropy, and can never be ne- 
glected. Several groupsSii^ii^ reexamined the Moriya' s 
theory and found expressions for the effective spin Hamil- 
tonian which includes both types of anisotropics. The 
linear SW theory applied to such models at T = al- 
lows one to obtain previously reported values of the spin- 
wave gaps at the centre of the 2D Brillouin zone, as 
well as to estimate the magnitudes of the anisotropic- 
exchange interactions. However, a detailed consideration 
of the model with the antisymmetric and symmetric DM 
anisotropics at nonzero temperatures is up to now absent 
from the literature. 

A very rough and simple approximation which can be 
used to study the effective magnetic model at finite tem- 
peratures is the mean field approximation (MFA). The 



MFA ignores effects of fluctuations and correlations be- 
tween the spins, hence, it fails for T near and gives no 
short-range order above the transition temperature. At 
very low T the noninteracting SW theory is useful, and 
it gives a successful prediction of the energy of low-lying 
excited states, and correctly reproduces the dominant 
term in the low-T magnetization. But, it fails near the 
phase transition point. To analyze the high temperature 
behaviour the 1/T expansion method can be employed. 
But, since the La2Cu04 crystal ordering temperature is 
much smaller than the magnitude of the superexchange 
interaction (Tn « J), the high-temperature expansion 
(to the first few orders in J/T) is not able to discuss the 
temperature region of interest, that is T near the transi- 
tion temperature. 

In the present paper time we consider the 2D spin- 
i anisotropic quantum Heisenberg antiferromagnet over 
the entire temperature range including both the symmet- 
ric and anti-symmetric DM interactions. We employ the 
technique of double-time temperature-dependent Green's 
functions within the framework of the random-phase ap- 
proximation (RPA). The first time such a scheme was 
used was by Tyablikov^ and he applied this formalism 
to the Heisenberg ferromagnet (the RPA for magnetic 
models is often referred to as Tyablikov's decoupling ap- 
proximation). This work was generalized by Liuii to ob- 
tain the longitudinal correlation function, and this latter 
study is important in the development presented in our 
paper. The important feature of this technique is that it 
deals with the entire temperature region and is in a good 
agreement with the SW theory at low-T, as well as with 
l/T expansions at high-T. In this paper, within such a 
scheme, we find the transition temperature at which long- 
range order would be established for an isolated plane. 
We obtain the excitation spectrum, sublattice magnetiza- 
tion and susceptibility tensor as function of temperature 
and coupling constants. We also employ the MFA and 
SW theories to compare results of all of these approxima- 
tion schemes, and note the essential differences between 
them. 

Of course, many investigations of the 2D spin-i have 
been completed previous to this work. We have already 
mentioned the most popular and simple methods to study 
spin models, that is phenomenological Landau theory, 
linear SW theory, the MFA. and high-temperature ex- 
pansions. They yield an analytical description of a wide 
range of physical properties and are very useful for prac- 
tical purposes. At the same time the great progress in 
the understanding of the ground state, thermodynamic 
properties, and spin dynamics of the Heisenberg magnets 
was made with the use of the newer and more compli- 
cated analytical schemes. Arovas and Auerbacbi^ used 
a path-integral formulation of the MFA theory within 
the Schwinger-boson representation. This method corre- 
sponds to the large- A'^ limit of the generalized SU(A^) 
model; however, various difhculties with this method 
have been discussed in the literatureii^iSS TakahashiSi 
has formulated and successfully applied the so-called 
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modified SW theory to the Heisenberg model which re- 
produced the results of conventional SW theory and is 
closely related to the Scliwinger-boson theory. For the 
one dimensional chain, Takahashi's modified SW theory 
yields very good agreement with Bethe ansatz results, as 
well as for the 2D classical ferromagnet at low-T (in that 
it agrees with Monte Carlo results). A self-consistent SW 
theory that is based on the boson-pseudofermion repre- 
sentation, was developed to study thermodynamics of 2D 
systems, and was also applied to S > 1 systems with an 
Ising-anisotropy 2D magnetsi^^ An important feature of 
all these methods is that they can be used to describe 
both the ordered and disordered {i.e. the case of no long- 
range order) states. 

Other related work includes: (i) The fermion rep- 
resentation to perform a 1/N expansion was used by 
Affleck and Marston;^ large-S" 2D Heisenberg antifer- 
romagnet in the long- wavelength limit; and (ii) based 
on the diagrammatic method for the spin operators, 
the thermodynamics and the longitudinal spin dynam- 
ics of Heisenberg magnets were studiedi^^iS^ However, 
the most note-worthy success in the investigation of this 
system is the work i of Chakravarty et alm^, who used 
a renormalization-group approach to the quantum non- 
linear cr model, the latter of which describes the low-T be- 
haviour of the 2d Heisenberg AF in the long-wavelength 
limit. 

As will become apparent below, the formalism that we 
have chosen to implement is more appropriate for this 
problem than any of those listed above, or the theories 
listed above arc too complicated to invoke when one goes 
beyond the 2D spin-i Heisenberg AF and includes spin- 
orbit couplings. 

The above few paragraphs summarize theoretical ef- 
forts that were directed towards the understanding of 
the 2D S=l/2 square lattice AF. The application of these 
and related work to describe the magnetic properties of 
so-called single-layer cuprate superconductors, such as 
La2Cu04, has attracted the attention of many theorists, 
and fortunately an extensive review of this work, written 
by Johnston, already existsiSl In this reviewSl one can 
find the comparison of the temperature dependence of the 
magnetic susceptibility for an AF Heisenberg square lat- 
tice calculated by different analytical methods and quan- 
tum Monte Carlo calculations, and, apart from the (post- 
review) data given by Lavrov et. al&, the application of 
the analytical predictions together with the numerical re- 
sults show very good fitting to the experimental data for 
the different single-layer cuprate compounds. 

Our paper is organized as follows. In S^we present the 
model Hamiltonian that we will study, introduce a con- 
venient coordinate transformation with which it is simple 
to complete analytical calculations, and then derive the 
transformation that relates the static uniform suscepti- 
bility in both coordinate systems. In ijlHI we derive and 
describe the MFA results, and then in ijlVl we present our 
derivations from applying the Tyablikov/Liu approach to 
our model Hamiltonian. In fj^Jwe present a detailed ex- 



amination of numerical results that follow from our work, 
including a comparison of MFA, RPA and SW theories. 
Finally, in HVll wc summarize our paper including a brief 
discussion of the remainder of the work that we have 
completed on the full three-dimensional problem. 

II. MODEL AND DEFINITIONS: 

A. Model Hamiltonian and the initial 
representation 

We consider a model for the Cu spins that are 
present in the Cu02 planes of a La2Cu04 crystal in 
the low-temperature orthorhombic (LTO) phase and 
employ a square lattice with nearest-neighbour inter- 
actions described by the following effective magnetic 
Hamihonianti^ 

= j^s,.s,+5]D,,.(s,xs,)+^Sr'r.y-s,. (1) 

This Hamiltonian consists of the superexchange inter- 
action together with the antisymmetric Dzyaloshinskii- 

Moriya (DM) interaction (D term) and the symmetric 

< — > 

pseudodipolar interaction ( F term) . As was discussed in 
the introduction, the DM and pseudodipolar anisotropics 
arise as a result of the mixture of Hubbard-type interac- 
tion energies and spin-orbit coupling in the low symmetry 
crystal structure. 

For the LTO phase, we use anisotropic interactions 
given by of the following form 

Da6 = ^(-1,1,0), D,, = -^(-1,-1,0), (2) 

and 
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where the corresponding coordinates, in what we refer 
to as the "initial representation" in the LTO phase, are 
shown in Fig. Qla). Note that the DM vector given in 
Eq. 121) alternates in sign on successive bonds in the a — 6 
and in the a — c direction of the lattice, as is represented 
schematically by the double arrows in Fig.QJb). 

We mention that the symmetric tensor F has been 
obtained by several author32iiii^iiSiS& in different forms. 
We have chosen the general form of this tensor, from 
which other specialized choices can be extracted. For 
instance, the form of the symmetric tensor obtained by 
Koshibae, Ohta, and Maekawai4 can be recovered from 
this definition if F3 = F2 — Fi. 

In the LTO phase the classical ground state is deter- 
mined uniquely^^ii^ and below the Neel temperature the 
Cu spin structure shows long-range antiferromagnetic or- 
der with weak ferromagnctism {viz. all spins cant out of 
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FIG. 1: (a) Coordinates in the initial representation, (b) 
Thin arrows — the Cu spins, open arrows — the DM vectors. 



the plane). To be concrete, in the classical ground state 
the spins are canted from in-plane antiferromagnetic or- 
der by a small angle given by 



1 

= — tan 



d/^/2 



J + 5(ri + r3) 



(4) 



and each plane has a net ferromagnetic moment in the z 
direction perpendicular to the Cu02 planes (weak fcrro- 
magnetism) . 

In the simplified case of the zero pseudodipolar interac- 
tion ( r 0) it was foundiSiSS that the ground-state spin 
configuration exhibits the rotational symmetry about the 
DM vector which is the origin of the Goldstonc mode 
in the spin-wave spectrum. Since in this simplified case 
there is a continuous symmetry in the ground state, the 
thermal fluctuations destroy the long-range order for any 
T > 0, according to the Mermin and Wagner theorenip^ 
In the general case of the model Hamiltonian of Eq. 
the continuous symmetry no longer exists and the spin- 



wave spectrum is gapped in the long wavelength limit 
q = 0. Consequently, the effect of fluctuations is reduced. 
That is, the DM (D 7^ 0) together with pseudodipolar 

( r 7^ 0) interactions can give rise to long-range order 
for low (but nonzero) temperatures even for the purely 
two-dimensional case (Tn > 0), and the Mermin and 
Wagner theorem does not preclude the possibility of a 
nonzero sublattice magnetization for nonzero tempera- 
tures in this general case. (Note that his does not imply 
that the transition to 3d long-ranged magnetic order is 
not influenced by the inter-planar exchange coupling, but 
simply that this latter coupling is not, in general, neces- 
sary to achieve such order.) 



B. Characteristic representation 

In solving this system, it is more convenient (theoreti- 
cally) to transform from the initial representation, given 
above, to the characteristic representation (CR) in which 
the quantization axis (z) is in the direction of a classical 
moment characterizing the ground state. In the present 
case there are two such classical vectors in the direction 
of the canted moments (recall that we are considering 
only a single Cu02 plane). Therefore, we introduce two 
rotated coordinate systems, as shown in Fig. |2J Spin de- 
grees of freedom in the initial representation are denoted 
by {Si}, but in the characteristic representation we use 
{(Ti}. (We follow the notation that i-sites belong to sub- 
lattice 1, whereas j'-sites belong to sublattice 2.) For the 
sites of sublattice 1 we apply a transformation of the form 



sm f 
sin^ 



sin 9- 
sin 9- 



V2cos6' 
V2cos6' 
2 sin 6* 





(5) 



and for sublattice 2 



The model Hamiltonian of Eq. ^ in terms of the new 
operators a reads 




The quantization axes (z) of the new spin operators at 
and cFj coincide with the unit vectors in the direction of 
canted moments moments Fig. |5| 



{i,3)c,b 

+ J2 {A{<^t^J +^i<^t) +Bo+a+ 

(7) 
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sentation, 



FIG. 2: Numbered arrows represent the Cu spin structure in a 
Cu02 plane. Two sublattices 1 and 2 are introduced. For each 
sublattice the spin coordinate system within the characteristic 
representation (i.e. after the transformations given by @ 
and 0) is shown. The thin net is shown only to simplify the 
visualization of the spin structure. 



where we introduced the following definitions 



Ji = J + Fi 



J2 = ^(ri - F3) + J (dV2) + [J + i(Fi + TsW, 



^3 = -i(ri - Fa) + ^ (dV2) + + ^(ri + r3)]2, 



Ja = —Ft sin 6 H — cos 9, 

V2 



. J1 — J3 „ J4 , .J1+J3 

, B = hi . 

4 ' 2 4 



(8) 
(9) 



The subscripts {i,j)ab and {i,j)ac in the summations of 
Eq. ((TJ imply the nearest neighbours in the ab and ac 
directions, as shown in Fig.^b). 

The form of the Hamiltonian in the characteristic rep- 
resentation is similar to an XYZ model, but is clearly 
more complicated since terms of the form a~aj — erf 
are present, which thus imply terms like Sf . Thus, we 
can extract from our results, in this representation, the 
magnetic susceptibility of the XYZ calculated in both the 
mean- field and random phase approximations, by setting 
the imaginary part B = 0. We will consider numerical 
results for this simpler model in a future publication. 



Magnetic susceptibility in the initial and 
characteristic representations 



We consider the response of the system, described by 
the Hamiltonian H in either the initial (Eq. (^) or char- 
acteristic representation (Eq. iQ), to an externally ap- 
plied constant magnetic field h. It is convenient to con- 
sider the application of this field in one direction only, 
which we take to be the a direction of the initial repre- 



N 



(10) 



1=1 



where a = a; or y or z, it is to be noted that a is not 
summer over in Eq. (|10|l . and N is the number of the 
lattice sites. 

The statistical operator of the system is required to 
evaluate ensemble averages of relevant physical quanti- 
ties, notably correlators and thermal Green's functions, 
and can be written as 



-^"Tr exp j/i" l^^"^^^'^''] ' 



where Si{t) = e^'^ Sie~^'^ is the operator in the Heisen- 
berg representation for imaginary time argument r, and 
Tr is the time-ordering operator. The zero-field suscepti- 
bility describes the response of the system to such a field, 
and is defined to be 



_ d{M° 



N N 



1=1 = l •^'J 



where 



N 



(12) 



(13) 



with correlators such as {TrS^{T)S" {0)) taken with re- 
spect to the zero field Hamiltonian H. 

The square lattice is bipartite and can be divided into 
sublattices 1 and 2. Then, by using the definitions 



Xii = 



X22 



X12 



2 



N/2 N/2 13 

/ (r.5;(r)5;,(0))dr, 

N/2N/2 p 
,= 1 j = l -^0 



2 

N 



X21 



N/2 N/2 p 



(14) 
(15) 
(16) 
(17) 



where 



i,i' E sublattice 1, & sublattice 2 

we can express the quantity of interest, x", as 
„ 1 



X 



:{Xn+X22 + X?2+X2i}- 



(18) 



6 



Then, using symmetry equivalent this simplifies (see be- 
low) the calculation of the zero-field susceptibility in the 
initial representation to 



X 



(19) 



The simpler form of Eq. ITJ vs. Eq. ^ makes clear 
that it is desirable to perform calculations first using the 
characteristic representation, and to then transform back 
into the initial representation. To this end we require the 
relevant form of the susceptibility tensor in the charac- 
teristic representation. To begin, let us perform transfor- 
mations Si = AcTi, S2 = B(T2 (-4 = [Oaa'], B = [5qq/]) tO 

the characteristic representation, such that the analogue 
of Eq. (Unjl is 

N/2 

1=1 
N/2 

- ^(^ox^tJ + bo^yG^ + 6„,a|)/i^. (20) 

i=i 

Note that we have generalized the applied field to be hi 
for sublattice 1, and h2 for sublattice 2, and in general 
we will treat these as two independent applied fields. If 
we define the components of susceptibility in the charac- 
teristic representation as 



N/2N/2 p 



Xll 



N 



i=l i' = l 



EE/(^-<W'^"'(0))d^' (21) 



^12'° ^nT.Y. / (T.<(r)a/(0))dr, (22) 
i=i j=i-'° 



then the susceptibility given in Eq. H14|) {N.B. in the 
initial representation) can be written as 



Xll — ,r 



N ^ dh? 

i—l ^ 
N/2 



dhf 



2 CT"^ CT^ CJ^fT^ CT^CT"^ 

+ 0-ayXll + aayCLaxXll + 0,ayO-azXll 
"I" ^azXll + O'ctz'laxXll + O-azO^ayXll j 



(23) 



and, in the same way (see Eq. (|19|l 'l 



X12 — ^ax^CixX\2 ~^^CiX^ayX\2 ~^^ax^azXl2 
-\- ClaybayX'12'^ '^(^aybaxXl2" ^O.aybazXl2'^ 
+ aazbazXl2'^' +aazbaxXl2'^" +aazbayXl2''''- (24) 



The quantities x'^ that are introduced above in 
Eqs. H21I22() have the following interpretation. For in- 

stance, the component ^12 determine the response 
of the expectation value '2/N^^J'^{af) of the spins of 
sublattice 1 to the magnetic field applied to the spins 
sublattice 2 (no field applied to the spins of sublat- 
tice 1) in the a' direction. Indeed, the perturbation 
H' = H — h2 X^jlfi '^f formally leads to the response 



J 



2 9"^/' 2 

^ 1=1 



N/2 N/2 p 



EE / (rr<(r)a/(0))drEEx?2°^°'- 

z=l , = 1 -^0 



(25) 



Similarly, the response of the spins of sublattice 1 to the perturbation H' = H — hf J2i=i given by 



2 d 



Ndh 



, N/2 N/2 N/2 p 

- E(^") -nT.T. (T.<(r)<'(0))dr ^ xf^"^ ■ 

1 1=1 i=li' = l-'0 



(26) 



r 



So, by substituting the inverse to the CR transforma- characteristic representation (which can be derived an- 
tion, given by Eqs. (|5l(j|) . into Eqs. 123I24|I . and taking alytically), one obtains the desired transformation be- 
into account that x°'''°'"=x'^'''^^=x'^"'^°^=X'^"'^''=0 in the tween the two representations, namely 
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X = Xii + xh 



2IX1I +X12 +X1I +X12 Xll X12 Xll X12 



sin2(6l) 
2 

+ cos2(0)(x?;-^-X?2^^), 

COS2(0) 



IXll +X12 +X11 +X12 +X11 +X12 +X11 +X12 J 

')(x?;"^-x?2"^), 

iXii X12 +X11 X12 +X11 X12 +X11 X12 ) 



(27) 

(28) 
(29) 



r 



III. MEAN FIELD ANALYSIS 

In this section we develop the mean field approxima- 
tion (MFA) for the system defined by Eq. and obtain 
the behaviour of the magnetic susceptibility and a defin- 
ing equation for the order parameter as a function of 
temperature. In part we include this derivation to make 
evident how the formalism of ^ III CI is applied to extract 
the zero-field uniform magnetic susceptibility. However, 
and more importantly, we will show that when the cant- 
ing angle induced by the DM couplings is small, there 
are significant deviations from the mean-field results, viz. 
quantum fiuctuation effects are large. Thus, here we es- 
tablish the MFA susceptibility with which to make these 
comparisons. 

Within the MFA we focus on one of the spins and re- 
place its interaction with other spins by an effective field. 
To this end the following replacement is used: 



(30) 



where a and b can be equal to any of x,y,z. It is to 
be noted that it is more convenient to perform the MFA 
calculations starting from the model in the characteris- 
tic representation, and thus we consider Eq. lO and the 
analogue of the above equation for the a operators. 

First, we find the equation for the order parameter. 
The Hamiltonian Eq. (7J) within the MFA reads as 



MFA 



(31) 



and we find that the order parameter, to be denoted by 
77, is found from the solution of 



14. 



itanh|^ZJ2((T- 



(32) 



where J2 is given by Eq. (jS]), and Z is the coordination 
number. From this equation it is immediately seen that 
within the MFA the Ncel temperature at which 77 vanishes 



rpMFA 



^(ri-r3) + ^(dV2) + [J+i(ri+r3)]2 

(33) 



Now, we find the susceptibility of the system within 
the MFA below T^^^^. First, we apply a magnetic field 
in the z direction of the sublattice 1 



i-sites e 1 sublattice. 



(34) 



The Hamiltonian within the MFA can be written as 

(35) 

where means sum over all sites i which are nearest 

neighbours of site j. Then 



af) = itanhJ^(^J2(a|) 
= itanh|^5]j2(a.f)| 



hi 



(36) 



We write the mean value of operators in the form 

= (aJ)o + = (^3^)0 + Sa'^, (37) 

where (erf )o = (c2)o — is the expectation value of cr^ 
operator in the absence of the field, and the term Scr^ is 
the part of (cr^) induced by the applied field. Since the 
applied field /if as well as the terms involving 5a^ are 
small, we may expand Eq. H36|l in powers of these terms. 
Then, we find 



Xll 



X^i^ 



6af 



fsech^ 



{|2J2.} 



ZJ2(f)2sech4{fZJ2r?} 
M 'n=o"l_(^)2scch4 I I^Jary}' 



(38) 



Due to the complicated couplings found in Eq. ((JJ), 

the transverse components are much more involved to 
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calculate. Applying a field in the x direction to the spins 
of sublattice 1 we consider 

H' = H -hl'^af, i-sites e 1-sublattice, (39) 



and within the MFA we thus examine 

i 

-E(f52^J + f52^J + f52^|)- (40) 
3 

Similarly, by applying a field in the y direction to the 
spins of sublattice 1 we consider 

i 

-E('^2^l + '^2'^J + f52^|)' (41) 



The solution of this systems, Eq. (|42|) turns out to be 



where 



f5?=5]{-2A(a;)+29i?(aJ)}, 



fjf =^{-2A(af)+23i?(af)}, 

where 3i? denotes the imaginary part of B. Then, the 
system of equations determining the transverse compo- 
nents of susceptibility Eq. and Eq. within the 
MFA scheme is found to be 



2 ^" 




- 35X21 - 


1 

2Z 


1^X21 




-^Bxf^"\ 




2 


- ^X2l"' 


-3sx2^"^ 




J2 „y„-: 

YX21 


- Axir 






J2 c'^crV 

-yXii 


- Ax2r 


-5sx2^^ 




2 


= Axll^' 






J2 cjVaV 

-yXii 


- Ax"2r 


-3sx2r'^"- 


1 

2Z 


J2 crWo-f 

9 ^21 


- Axir 


-3i?x?^'^^ 





Xii — X22 ^ Xii — X22 



X12 — X2I — X12 ~ X2I 



Xii — X22 — Xii — X22 



X12 — X21 — X12 — X21 



J2/2+A 


J2/2-A 




4Zw2 


J2/2+A 


J2/2-A 






35 / 1 








35 / 1 






^5/ 



(43) 



where 



Loi = y^(J2/2 + A)2-352, 



UJ2 - V (-^2/2-^)2 -35^ 



(44) 



Using the relation between the components of susceptibil- 
ity in the initial and characteristic representations given 
in Eqs. I|27() - (|29() . we obtain the final result for zero-field 
uniform susceptibility within the MFA below the MFA 
ordering temperature, 

1 1 



MFA 



X 



MFA 



AJ1 + J2 



(45) 



y MFA _ 



1 sin2(6') cos2(6') 



sech" 



{IZJ2V} 



4 J2 - J3 



MFA 



1 cos2(6i) sin^(6l) 



T+J2 sech^ {fzJar?!' 

(46) 

sech^ ||zJ2?7| 



(42) 



(47) 

with the equation for the order parameter ry given by 
Eq. (For d = T,^0, implying that 9 = and J2 = 

J, the above seemingly complicated results indeed reduce 
to the correct MFA expression for the susceptibility.) 

The following comments on the MFA result are in or- 
der. First, note that for physical values of d and F^ 
(d, Ti <C J) the canting angle out of the xy plane 
is very small; thus, since the AF moment is in the yz 
plane and nearly aligned along the ±y axes, x^ diverges 
at T^^^, but the other two components remain finite at 
the transition. However, while the x component of the 
susceptibility remains independent of temperature, since 
the canting produces a net FM moment in the z direction 
that is coupled to the y component of the local moment, 
there is an additional increase of x^ as the transition is 
approached from below. 

Now consider the paramagnetic temperature region 
(T > Tn), for which the only components with nonzero 
spin expectation values are those driven by the applied 
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field. Following similar considerations to above, the final 
results for the components of susceptibility in the initial 
representation for high temperatures (T > Tn) reads 



r 



MFA 



1 1 



4,/i + r' 

2 



X 



X 



y MFA 



MFA 



lsin^(6') lcos2(6l) 
4T- Ja ^ 
1 cos2(6') 



4 T + J2 
1 sm^{e) 



AT + J3 AT-J2 



(48) 
(49) 
(50) 



Note that in the hmit T ^ T^'^^ = J2 we obtain that 
the x,y components of the susceptibility are continuous 
at the transition, whereas the z component of the suscep- 
tibility diverges at the Neel point, from above or below, 
owing to the presence of the weak ferromagnetic moment 
that first develops at the transition. 



IV. LINEAR RESPONSE THEORY WITHIN 
THE RPA 

A. Susceptibility below Tjv 

In this section we derive expressions for the static, uni- 
form susceptibility within the RPA below the ordering 
temperature, Tpf. Note that this temperature is deter- 
mined with the RPA, and is not equivalent to that found 
in the previous section. 

We employ thermal Green's functions in the analysis 
of the spin Hamiltonian given in Eq. (^) with spin i. The 
definition of such Green's functions for two Bose opera- 
tors A, B and the corresponding equation of motion, arc 
given by 



Gab{t) = {TrA{T)B{0)) 

dGAB(T 



(51) 



dr 



5{T){[A,B]) + {T^[H{r),A{r)]Bm. (52) 



As discussed in the introduction, we adopt a procedure 
that was introduced by Liu^ as this technique allows for 
us to find longitudinal component of the susceptibility. 
To this end, we introduce the perturbed Hamiltonian (in 
the characteristic representation) 

i 

where / is a small fictitious field; note that the field is 



applied to the spins of suhlattice 1 only, and within the 
present paper we restrict / to be constant and static. 

In the imaginary-time formalism, the Green's functions 
to be used are 

GLM = (rrT+(r)<7-(0))/, 
GI~{t) = {TrCr^{T)a-{0))f , le sublattice 1, 
Gk(r) = (r.a+(rX(0))/, 
Gi'niT) = {Tr(T-,{T)a-{Q))f , sublatticc 2, 

(54) 

where the expectation values are taken with respect to 
the perturbed Hamiltonian in Eq. (|53|l . After an expan- 
sion in a power scries of / we can write 

Gfjr) = Gf^{r) + fG[i\r) + 0{f). (55) 

Since g[^{t) = Ginir), from now drop the superscript 
and use 

Glir) = G,„(r) + fGg\r) + 0{f). (56) 

Also, we introduce 

K(t))^ = «)+./v, + 0(/2), (57) 

where, due to the translation periodicity (erf) = 77, the 
order parameter at / = 0. 

The equation of motion for the Green's function G/„(r) 
is given by 



dr 



+ {Tr[HcR{T),a+{T)]a~{0))f - /Gf„. (58) 



In order to solve this equation for the Green's function 
it must be linearized. We will use the random phase 
approximation (RPA) , in which the fluctuations of cr^ are 
ignored and the operator is replaced by its mean value 



this is the so-called Tyablikov's decoupling. 
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For example 

(T.af(r)a+(r)a-(0))^^ 
^ (af(r))/(r.a+(r)a-(0))/=(af(r))/Gf^.(r). (59) 

After this decoupHng is introduced, Eq. lf5S|l is found to 
be 



dG^ (t) r ~i 

= 2<5(r)<5,„(af )/ - J2 W)'' [^Gf,^.)„(^)-SG;,;,)„(r)] + M<yUsir)y Glir)} 

Sab 

-5] {2(af(r))/[AG;^,)„(r)+i?*Gj;,)„(r)] + M^Usir))' G{jr)} JGUt), 



(60) 
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where ^ refers to a summation over the nearest neigh- 
bours of the site I in the ab direction, and similarly for 
— see Fig. mb). Here, all sites I + S belong to the 
sublattice 2. 

We introduce the Fourier transformation in the 
momentum-frequency representation for the Green's 
function and the spin operator 



(61) 



h.rn 



-ifc Ri 



(62) 



where the expansion in Eq. 1)57(1 and the linear response 
to the uniform perturbation expressed by vi(fc) — (5(fe)vi 
were taken into account. In the transformation given by 
Eqs. H61I62|I . the sum over fe runs over points of the 
first zone in the momentum space, and = 27rn//3 for 
n £ Z are the Bose Matsubara frequencies. Then, we can 
write down the equation for the Green's function (t) 
in the form 



- /G(2(fc,Wm) 

- 2J2[?? + ./v2]G(2(fc,W™) 

+ 2ZSfc[77 + /vi]G4-(fc,c^„0, (63) 



where, as before, Z is the coordination number, and we 
introduce 



(64) 



Ik 



-{coskx + cos fey), 7^ ~ —{coskx — cosfcy). 



From these we can write down the following two equa- 
tions: 



-G 
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-Gi2 + AkG22 — -BfcG22, 



(65) 



i|^„(i) 



1 



-Gi2 



2Z7] 

V2 0/2. 

77 2 



- — -z-^^12 + -7r^l2 



■h. 

2 



Vl 



V 



—AkG22 + AkG^^^ 
(66) 



where in all equations we drop the wave vector and fre- 
quency dependencies for the Green's functions, that is 
G = G(fc,tJ,„) and G^^) = G^^) (fe, tj,„). 

In the same way we obtain the equations of motion 
for the other Green's functions (see Eq. (|54|l ) within the 
RPA scheme. The final systems of equations for zcrotli- 
and first-order quantities can be written as 



/ i^m _ J2 
\2Zr] 2 



G12 — ^fcG22— SfcG22, 



iuJm, J 2 



2Zt] 



-Tr ] G22 — AkGi2 — BkGi2 — ^5: 



ia7„; J 2 

2z7] ^ y 

1 / iuj„i , J2 



G12 — -AkG22 + BIG22, 



\2Zr1 ^ Y ) ^22 - -AkG^2 + BIG12] 



(67) 



- 1) - ^^^^^ - + {?f + - f ) ^ 2^;}°- 



(68) 
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JA ^(1) 

JA ^(1) 



(1)- 



ILOm 



f ) G[r 



(1)- 



V2 ^ Vi / i^Tn _ 

r; 2 77 V2Z77 

vi J2 _ / i^m 
?/ 2 rj \2Zri 

V2 ^ _ vi / i^m 
?7 2 ?7 \2Zr] 



2 
2 



'12, 

2Z^ 

G121 
1 



G 



2Zri 



22j 



g; 



(69) 



where we have taken into account the relations 

G12 = G21, Gil = G22, 

G\2 ~ G21, Gil '-'22- (70) 

The poles of the zero-order Green's functions G have 
to be the same as the poles found for the first-order ones 
G^^-' . This can be seen directly by comparing the struc- 
ture of the systems of equations for the corresponding 
quantities: the system in Eq. (|67|) for the zero-order func- 
tions is identical with the systems in Eqs. H68I69|1 for the 
first-order ones, except for the free terms. The free terms 
in the first-order systems are determined by the zero- 
order Green's functions, thus, the first-order quantities 
can be written down in terms of the solution for the 
zero-order system of Eq. (|67|) . and the as yet unknown 
quantities Vi and V2. 

To calculate vi_2 we use a relation connecting v and 
the Green's functions G^^^k,0~). From the definitions 
in Eq. (|54|l and the expansion in Eq. H57|) we have 

Gf,(0-) = i-«)^ = i-r;-/v,;, (71) 



while the expansion in Eq. lf5B|l leads to 



GiiQ-) = G,,(0-) + /G«(0-) = 2 - ^ + /Gil'(O-). 

(72) 

Thus, we can write down — v,; = g[1\o^), and after 
Fourier summation one obtains 

-^1 = |E^n(fe,0-), (73) 

k 

= ^J2^^^2ik,0-). (74) 

k 

The solution of the system in Eq. H68(l gives us the 

first-order Green's function G22'' [k, lu„i) and therefore V2. 
Similarly, to find vi we use Eq. (|69ll . 

The solution of the system of equations in Eq. (|67|l for 
the zeroth-order Green's functions turns out to be 



Gi2(fe,LJ„) = 
G22(fe, W,i) = 

Gi2{k,UJn) = 
G22(fc,W„) = 



J2/2+Ak 



iOlik) J 


\un-e\{k) 


•hn—Ak 




LJ2(fc) , 


/ \UJn-e2{k) 


•hn-^Ak\ 


1 


iOlik) J 


it^„-ei(fc) 


.hl^-Ak 


\ 1 



-iBl 



LJ2(fc) / \U:n-S2[k) 
1 



1 



1 - 
1 - 

1 - 
1 - 



J2/2+Afc^ 



c^l(fc) ) 




fei(A;) 


J2l2-Ak 




1 


L02{k) 


) it^7 


i+£2(fc) 


J2l2+Ak\ 




1 


Wl(fc) ) 




f£l(fc) 


•hl^-Ak 




1 


LU2{k) 


) it^r 


i+£2(fe) 






1 


+ LJ2ik) ( 


iuJn- 


-£2(fc) 






1 


^2(fc) ( 




-£2(fe) 



(75) 



12 



where the spectra for the out-of-plane £i(fc) and m-planc e2{k) modesi^ are given by 

e2{k) = 2Zriuj2{k) = 2Zri^ {J2/2 - A^f - \Bk\'' . 

(76) 

After the substitution of the results in Eq. H75|l into the system of equations in Eqs. H68I69|I . and then using the 
solutions for G^^^'' (fc, Wm), (fc, lj„i), the results for quantities vi and V2 are found to be 

1 + 477^7261 

V1 + V2 = - — ^-^7=r. (78) 
1 - 877-^ Ca 



where 



r 2^^f/ (^j^/2)^.Al-\B^\' \ nie,)-n{s2) ( {J2l2Y-Al^W \ n{eMn{e2)+l \ 
' ^fc^^V ^i(fc)^2(fc) y £2(fc)-ei(fc) V ^i(fc)c^2(fe) y ei(fc)+e2(fc) / 

^ 2 ^ / (J2/2+^fc)2/3/2 |gfcH2n(£i)+l] (72/2-^^)^/3/2 \Bk\''[2n{e2)+l] \ 
' A^4'\c.2(fe)sinh2^ c.2(fc)£i(fc) +c.|(fc)sinh2^ ^Uk)^2{k) j' 

2 / (J2/2+ylfc)/3/2 ^ (J2/2-Afc)/3/2 | r ^ r ^^^^ n-i 

= ^E|^-JT^ + ^-JT^|- heren(e,2) = [exp(/?s.2(fe))-l] . 



(79) 



Now let us find the quantities which determine a linear 
response to a magnetic field applied to the one of sublat- 
tice - see Eqs. (|25I26II . The longitudinal z components 
of the susceptibility in the characteristic representation 
are given by 



All — 



Of 



f=o 



= vi, 



A12 ^ 



Of 



- V2, 

(80) 

where the expansion of Eq. (|57|l was used. The trans- 
verse X and y components of the susceptibility tensor are 
determined in the terms of Green's functions as 



All 



N 



^12''"' = |E /Vr^r(T)<'(0))dr, 



(81) 



where a = x,y. By substituting the solutions in Eq. H75() 
into the definition in Eq. (|81|1 for the transverse compo- 
nents of susceptibility, we easily obtain exactly the same 
result that we have already found within our MFA cal- 
culations - that is, Eq. (|^ . 

Then, using Eqs. H27(l - (|29|l the components of the sus- 
ceptibility in the initial coordinate system of Eq. are 



found to be 



1 1 



4 Ji + J2 ' 
1 sin2(6l) 

4 J2 - J3 
1 cos2(6') 

4 J2TJ3 



-Kcos2(6')[vi - V2] 



-sin^(6')[vi +V2] 



(82) 
(83) 
(84) 



For completeness, we mention that we have also per- 
formed the theoretical investigation of this model 
within spin-wave (SW) theory, and the final result for 
the components of static susceptibility turns out to be 



X 



X SW 



1 



,V SW 



4.J1 + J2 

1 sin2(6i) 



X 



z SW 



1 COS^(e) , ,,;„2/z,Nc2. 



4 J2 - J3 

d2 



4 J2 + Js 



+ sin^{e)S^C2 



(85) 
(86) 
(87) 



It can be noted that the difference in the results within 
the RPA, Eqs. |jHlll-l|Hl|), and spin- wave theory, Eqs. 
(|87|1 . came from the calculation of the components of the 
susceptibility in the direction of the sublattice magne- 
tization (that is Xii'^ and Xu"^')- The spin- wave the- 
ory gives unity in the denominator of the expressions 
for xfi'^' and Xi2'^" in Eq. and 5=1/2 instead 

of the order parameter rj everywhere in the numerator. 
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The similar situation takes place for antiferromagnetic 
Heisenberg model within the RPA schemeii and spin- 
wave theoryi^ 

We also mention that the transverse components of the 
susceptibility in the characteristic representation 143|) are 
equal within the MFA, RPA, and SW theories. 



B. Related Thermodynamic quantities 



C. Susceptibility in the paramagnetic case 

When the temperature of the system is above the Neel 
temperature, T^v, there still exists short-range magnetic 
order. To model such an order— we introduce a ficti- 
tious field h pointing in the direction of the sublattice 
magnetization, that is the z direction in the characteris- 
tic representation. To this end, the Hamiltonian 



In order for the above RPA theory to be complete, we 
need to determine the behaviour of the order parameter 
and the transition temperature. 

The above expressions for the components of sus- 
ceptibility Eqs. H83I84|I . and for the elementary excita- 
tions (spin waves) given by Eq. (|76|) . include the as-yct- 
unknown value of the order parameter 77. From the defi- 
nition on the Green's functions wc can obtain 



Gnn{T = ) = (cr„cr+) 

2 



(88) 



where G„„(0-) - ^ G22(fc, 0") 



Substituting G22(fc, w) from Eq. (|75|1 . and performing the 
summation on the Matsubara frequencies, the equation 
on the order parameter turns out to be 



N 



E 



J2/2+Ak 

LJl(fe) 

, J2/2-Ak 



UJ2{k) 



[2n(£i)+l] 
[2n{e2)+l] 



(89) 



Since order parameter (|89|1 (sublattice magnetization) is 
temperature dependent, it follows that the spectrum of 
elementary excitations (Eq. H76|) ) is also temperature de- 
pendent. 

The Neel temperature at which 77 vanishes within the 
adopted RPA approximation is determined by 



J2/2+Ak , Ji/'^-Ak 



(90) 



By putting 7] we can find that z-componcnt of sus- 
ceptibility in Eq. (jHH) 



X 



.-sin^(.)-^^^l-° 



4 J2 + J3 



I-TTn ' 



(91) 



diverges at the Neel temperature, whereas other compo- 
nents of susceptibility remain finite as the Neel point is 
approached from below. 



(92) 



is used, and the limit h is taken after the calculation 
is carried out. To obtain the susceptibility above the 
Neel temperature, it is convenient to introduce an order 
parameter defined by 



y = lim {2Zr]/h). 



(93) 



The calculations for the model are very similar to the 
ones above presented. It is easy to show that paramag- 
netic version of the equation on the order parameter in 
Eq. ^ leads to 



E 



1 + y( J2/2+Afc) 



N^Zp{ (l+2/(J2/2+Afc))2 



y{J2/2-Ak) 



il + y{J2/2-Ak)r-y^\Bk\ 



(94) 



The quantity y approaches to infinity as the temperature 
is lowered to Tn- Indeed, putting y — > 00 in Eq. H94() 
we find the temperature at which y diverges, which is 
nothing but Neel temperature. 

By a procedure similar to the above presented (that 
is, the RPA scheme below Tn) the different components 
of the magnetic susceptibility in the paramagnetic phase 
are found to be 



1 



1 

4 Ji + J2 + 2/y' 
1 sm^{9) 
4J2-J3 + 2/y 

1 cos2(6') 
4 J2 + J3 + 2/y 



(95) 



COS^{0) 



+ sin2(6i) 



1 + 8y{l+yJ2/2)Di 

(96) 

2 



y'D2 



1 - Sy^Ds 



(97) 



where 
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1 2^ {i+y{^j^l2+Au)){l+y{J2/2-Ak))-y^\Bk? 
2Z2/3 TV ^ {(l+y(J2/2+Afc))2 - y^\Bk?}{{l+y{J2l2-Ak)f ~ y^\Bu?} ' 

{l+y{J2/2+Ak)r + y^\Bk\' il+yiJ2/2-Ak)r + y^W ' 

{(l+2/(J2/2+^fe))2 - y^Wr ^ {{l+yiJ2/2-Ak)r ~ yWV\ ' 

l+y{J2/2+Ak) l+y{J2/2^Ak) 
(l+y(J,/2+Ak)r-y^\Bk\^ + (l+y{J^/2-AkW-y^\Bk\\ ' 

(98) 



Do 



1 



2Z^PN 



E 



1 



2Z'^p N 



E 



By putting ?/ — > oo we obtain that the components of sus- 
ccptibihty ^-nd ^^'^ continuous at the Neel point, 
whereas the z-component of susceptibihty diverges in the 
y ^ oo limit at the Neel point, the latter result reflect- 
ing the presence of the spontaneous canted ferromagnetic 
moment in the z direction. 



D. Susceptibility in the T — limit 

As we will present in the results discussion, the dimen- 
sionality of the parameter space that seems to be relevant 
to the cuprates is large, but there are only a few impor- 
tant values that determine the physical properties of the 
system. Here we discuss two key experimentally obtain- 
able quantities, and their relation to the above theory. 

It has been reported, using inelastic neutron scattering, 
that the out-of-plane (ei) and in-plane (22) spin- wave 
gaps are 5.0 and 2.3 meV, respectively, in the LTO phase 
of La2Cu04 crystalu^ Using these results let us predict 
the ratio of the components of susceptibility x^ /x^- The 
zone-centre (fc = 0) spin-wave gaps are given by 



£1 - 2?7V(J2 + Ji)(J2- J3), 

62 = Zt]^/{J2- Jl)iJ2 + Js), 



(99) 



and they arc real if J2 < Ji, Ja- So, from these relations 
we obtain 



£2 < £1 ^ -^1 > J3 



(100) 



Also, in the T = limit the y component of the suscep- 
tibility in Eq. (|46|l is given by 



then 



y MFA ^ l_sin^(^ ^ i A_A 
4 J2 - J3 4 J| - J| 

y MFA j2 _ j2 



MFA 72 _ 72 ■ 
A '^2 "^3 



Therefore, within the MFA 

X MFA ^ y MFA ^ 7^7 
X >X ^ Jl > J3. 



(101) 



(102) 



Thus, if S2 < £1 (£2 > £1); in the limit of zero tempera- 
ture the MFA predicts that x^ < ix^ > x"")- 



In the limit of the small anisotropy d, F <g; J the com- 
ponents of the susceptibility at T = within the MFA 
turn out to be 



X MFA ^ z MFA 
X ^ X 



1 

8J' 



X 



y MFA 



32J2(Fi-F3)' 
while the expressions for the spin-wave gaps are 



ei « Z?/v/2J(Fi -F3), 62 « Zi]d/V2. 



(103) 



(104) 



We can see that components x^'^ ^re almost independent 
of the anisotropy parameters, while the x^ component is 
very sensitive to the ratio between the anti-symmetric d 
and symmetric Fi — F3 parameters of anisotropy. Then, 
the ratio between the components of the susceptibility is 
given by 



^,x,z MFA I _ 
X I 61 



X 



y MFA 



62 



(105) 



It can be noted that within the MFA scheme the dif- 
ferent components of the susceptibility, i.e. x^? X^ ^-^d 
X^, are determined by the contributions from the trans- 
verse components of the susceptibility in the character- 
istic representation. Indeed, as should be expected, the 
longitudinal components of the susceptibility in the CR 
(see Eq. (jSHJ) are equal to zero in the T ~ limit. As 
showed earlier in this paper, in the characteristic rep- 
resentation the RPA and SW theories lead to the same 
result for the transverse components of the susceptibil- 
ity as the MFA does. Since the longitudinal components 
in the CR, given by the Eq. within the RPA, and 
their simplified expressions within the SW theory (see 
Eqs. H86l87|l ') become negligibly small in the T = limit, 
wc predict that RPA, SW and MFA within the reason- 
able range of the model parameters {d, F ^ J) satisfy 
the ratio of Eq. H105|l . and the different components of 
the susceptibility at T = can be approximated by the 
Eq. fTU^ . 

We also note the analogy with the pure 3D Heisenberg 
model where, in the limit of zero temperature, all ap- 
proximations considered here give the same magnitude 
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for the transverse components of the susceptibihty, and 
zero for the longitudinal one^ 

V. RESULTS OF CALCULATIONS 

In this section we present the results of a numerical 
investigation of the magnetic properties of the system 
modelled by the Hamiltonian given by Eq. ^ based on 
the above presented analytical formulae. Specifically, we 
arc interested in the temperature dependencies of the var- 
ious components of the susceptibility for different values 
(specifically, ratios) of the model parameters. Further, 
we will numerically demonstrate the correlation between 
the magnitudes of the two spin-wave gaps in the exci- 
tation spectrum and the behaviour of the susceptibility 
components; this relation was discussed analytically in 
the previous subsection. Also, and most importantly, we 
will make clear the role played by quantum fluctuations 
by comparing the results of the different approximation 
schemes. 

In what follows, we will mainly examine one set of pa- 
rameters that are suggested from experimental measure- 
ments discussed in the previous subsection, and this will 
allow us to "zero in" on a parameter regime. However, 
since have developed the theory for one plane and not 
a 3d solid, we do not necessarily expect this set of pa- 
rameters to be representative of a system like La2Cu04; 
instead, as we discussed in the introduction to this pa- 
per, this approach will allow us to determine if a one- 
plane approach is adequate, since, as we and others have 
discussed, a true Tc > phase transition is possible for 
one plane and thus could possibly be sufficient for this 
system. 

In the present calculations wc express all model param- 
eters in terms of J. Also, as will be made clear below, 
instead of using the set of parameters Fi, F2, and F3, we 
deal with combination of parameters Fi — F3, Fi -I- F3, 
and The chosen magnitudes of the model parame- 
ters d and F give the reported magnitude of gaps in the 
spectrur 
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£0 = ^1 ~ 5 meV, Ej = £2 ~ 2.3 mcV, (106) 

at the temperature T = rjv/3 for the superexchange 
value J = 130 meVi^ This leads to the parameters given 
by 

d/J = 0.02, (Fi - F3)/J = 0.42 x 10"^ (107) 
where, as discussed below, we have set 

(Fi -I- F3)/J = 0, F2 = 0. (108) 



A. The AFM order parameter, spin-wave 
excitations, and Tjv 

To examine the different analytical schemes used in the 
previous sections, we compare the representative solu- 
tions of the order parameter, 77, within the RPA method. 



Eq. and within the MFA, Eq. Most interest- 

ingly, we note the results shown in Fig. I^a) look very 
similar to the corresponding ones for the pure 3D quan- 
tum Heisenberg antiferromagnet within the RPA and 
MFA schemes 

Since the order parameter ry is temperature dependent, 
it follows that within the RPA scheme the spin-wave 
spectrum (sec Eq. (|76f) ') is also temperature dependent. 
In Fig. Ofb) we present the behaviour of both modes 
in the excitation spectrum at the long wavelength limit 
(fc = 0) (energy gaps) with respect to the relative tem- 
perature {T/Tn); these results compare favourably with 
the experimental measurements of the same quantity^ 




(b) 

FIG. 3: (Color online) The (a) order parameter vs. T/Tn 
within the RPA method (black solid line) and the MFA (red 
dashed line), and the (b) spin- wave gaps, in units of J, in the 
spectrum of elementary excitations vs. T/Tm within the RPA 
method. In both of these figures we have used d/J — 0.02, 
(Fi - Fa)/ J = 0.42 x 10-^ Fi -f Fa = 0, and Fa = 0. 

Now let us show that in contrast to the MFA ap- 
proach, where Tn = J2 ~ J is almost independent of 
the anisotropy, the Neel temperature within the RPA 
analytical scheme is very sensitive to model parameters 
d and Fi — F3. Figure^ shows the zero-temperature en- 
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FIG. 4: The (a) T = energy gaps, in units of J, vs. the 
DM parameter d/J, as well as (b) the Neel temperature Tjv, 
in units of J, vs d/J. In both of these figures, (Fi — Vs)/ J = 
0.42 X 10"^ Fi + F3 = 0, and Fa = 0. 



ergy gaps and the Neel temperature as functions of the 
DM antisymmetric exchange interaction d/J within the 
RPA method. As one can see, the energy gap ei is almost 
independent of the d/J, while £2 depends almost linearly 
on the DM interaction d/J, and in fact goes to the zero 
in the limit d/J ^ 0. As a result, when d/J = the 
Goldstone mode appears in the spin-wave spectrum and 
thermal fluctuations destroy the long-range ordering for 
any T > 0. Consequently, the Neel temperature drops to 
zero in case of d = 0. 

In the next two figures we present the dependencies 
on the model parameter (Fi — Fa)/ J. Now, the en- 
ergy gap £2 is independent of the parameters of sym- 
metric anisotropy F and thus determined by the DM 
interaction d/J alone, while the gap £1 varies strongly 
with (Fi — F3)/J. As in the above case, in the limit 
of Fi — F3 — !■ the mode £2 in the spectrum becomes 
gapless and, therefore, the transition temperature to the 
long-range ordered state would be suppressed to zero. 

One can understand the above results for the zone- 
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FIG. 5: The (a) energy gaps, in units of J, vs. (Fi — F3)/J, 
and (b) the Neel temperature, in units of J, as function of 
(Fi-F3)/J. 



centre excitation spectrum immediately from the the ex- 
pressions Eq. (|104|1 in the limit of small anisotropy. Our 
numerical results, shown in the previous figures, demon- 
strate that these expressions are valid over a large range 
of parameter values. 



B. Parameters regimes 

We now summarize our numerical results with regards 
to the dependence of various thermodynamic quantities 
on the material parameters appearing in the Hamilto- 
nian. 

Firstly, we find that the Neel temperature is almost in- 
dependent of the (Fi -I-F3)/ J within the reasonable range 
of the model parameters (see below). In fact, in order to 
argue for the independence of the thermodynamic quan- 
tities central to this study on certain material parameters 
that appear in the Hamiltonian, viz. F2 and Fi -I- F3, in 
Fig. El we show two representative plots for the order pa- 
rameter and the susceptibility within the RPA scheme 
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FIG. 6: (Color online) The (a) order parameter vs. T/J 
for the different values of r2/J and (Fi + Fs)/ J, and (b) the 
susceptibility, in units of 1/J, vs. T/J for the different values 
of r2 and Fi +F3 discussed in the test. In these plots we have 
fixed d/J = 0.02, and AT/ J = (Fi - rs)/J = 0.42 x 10"^). 



for the constant values of the d and AF = Fi — Fa (again 
in units of J). That is, in each of the plots in Fig. |H|we 
have simultaneously plotted ten data sets each with the 
different values of the F2 and Fi +F3, where the parame- 
ter ratio F2 has been varied from the value —10^ x AF up 
to the 10^ X AF, and Fi + F3 from the value -10^ x AF 
to the 10^ X AF (all in units of J). As one can see, even 
for such a large range of the parameters, one can hardly 
see the difference of the absolute values of the Neel tem- 
perature, order parameter and susceptibility. 

Thus, to study the magnetic properties of the system 
we can use only the DM interaction d and the combina- 
tion Fi — F3 of the symmetric tensor components as two 
independent parameters, and so wc conclude (similar to 
others^SiS) that the system can be studied using F2 = 
and Fi -h Fa = 0. 

In various limits, it can be shown that this result fol- 
lows from the above presented analytical work. The non- 
diagonal term F2 of the symmetry anisotropy tensor is 



involved in all expression through the combination in the 
J4 (Eq. (jSJ). For the reasonable anisotropy parameters 
(that is F < d <C J) the spins are canted by a very 
small angle 9 ^ d/J, and as a result we can neglect the 
term F2 sin 9 ~ dT2 / J with respect to d, and hence we 
can ignore the quantity r2 in all our formulae. Simi- 
larly, Fi + F3 is involved in the formulae through the 
combination in the J2, J3, and canted angle 9, where it 
appeared only as the combination J+ i(Fi +r3). Thus, 
the parameter Fi -I-F3 can be ignored with respect to the 
superexchange interaction J (sec Eq. (jHJ). 

Therefore, one can assert that the model Hamiltonian 
of Eq. ^ leads to the same results as, for instance, the 
model described by the spin Hamiltonian 

H = ^[JS, • S, - AF^fS-; + D„(S, X S,)], (109) 
where we define AF = Fi — F3. 



C. Susceptibility 

Now let us consider the main focus of our paper, that 
being the behaviour of the different components of static 
uniform magnetic susceptibility as a function of temper- 
ature. Our results for x^, for the parameters discussed 
in the previous subsections, are shown in Fig. [7| (recall 
our earlier result that the MFA and SW theories predict 
the same T- independent value for this quantity). The x 
component of susceptibility below the Neel temperature 
is the temperature independent and is equal to w 1/(8 J) 
within the MFA (Eq. (gSJ), the RPA scheme (Eq. l|g^ ). 
or spin-wave theory (Eq. (|85|l ). However, above the or- 
dering temperature, the RPA and MFA yield different re- 
sults, with a weak T-dependence within the RPA, while 
a strong Curie-like falloff is found within the MFA. 
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FIG. 7: (Color online) The susceptibility x^i in units of 
1/J, within the RPA (black solid fine) and MFA (red dashed 
line), for the parameters values d/J — 0.02 and AF/J = 
0.42 X 10~^. Below Tjv these theories both predict the same 
constant value that is independent of temperature. 



As we will discuss in a future publication, this be- 
haviour changes if one includes 4-spin ring exchange, or 
goes beyond the Tyablikov RPA decoupling scheme that 
we employ in this paper. This is important since the ex- 
perimental data of Lavrov et al.^ shows a small nonzero 
slope of vs. T. Indeed, a successful comparison with 
the small slope seen below Tjv in experimental data^ nec- 
essarily requires that we go beyond the treatment of spin- 
wave interactions and/or Hamiltonian that are included 
in this paper. We emphasize that the necessity of go- 
ing beyond the Tyablikov RPA decoupling to obtain this 
slope is a manifestation of the presence of strong quan- 
tum fluctuations, a theme that will be repeated in our 
discussion in this and the next section of this paper. 

Our results for the y component of the susceptibil- 
ity, x^, are shown in Fig. |S1 These plots show that be- 
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FIG. 8: (Color online) The (a) susceptibility within the 
RPA (black solid line) and MFA (red dashed line), and (b) a 
comparison of the RPA (black solid line) and spin-wave (SW) 
(blue dotted line) results below Tjv. As in previous figures, 
we are using d/J = 0.02 and AP/J = 0.42 x 10"^. 

low the ordering temperature the RPA scheme leads to 
the good agreement with the MFA scheme near the T/v 
(O.ST/v < T < Tat), and good agreement with the SW 
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theory at low-T (that is, for T < Tn/2). Above the Neel 
temperature, the RPA and MFA theories lead to very dif- 
ferent results. The MFA method gives an abrupt decrease 
of to a value that is close to that of the purely trans- 
verse component x^ ~ l/(8>/) (see inset of this figure), 
while the RPA leads to a much more gradual decrease of 
the value of x^ with the temperature. 

The z component of the susceptibility, x^ i is shown in 
Fig. ini We find that at low T, as was also found for x^, 
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FIG. 9: (Color online) The susceptibility within (a) the 
RPA (black solid line) and MFA (red dashed line), as well as 
(b) a comparison of x^ below the ordering temperature within 
the RPA (black solid line) and spin-wave (SW) (blue dotted 
line) theory. As in previous plots, we have used d/J — 0.02 
and AP/J = 0.42 x 10"^ 

the RPA is in the the good agreement with spin-wave 
predictions. Near the transition temperature, the RPA 
method leads to the qualitatively different behaviour of 
the z component of susceptibility with respect to both 
the MFA and spin- wave formalisms. 

The differences between the MFA vs. RPA data shown 
above can be understood using the following reasoning. 
Firstly, consider below the Neel temperature. The canted 
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moments which develop are confined to lie in the y — z 
plane; as well, they are ferromagnetically ordered in the 
z-direction (recall that we are studying a single Cu02 
plane). Then, one can see that within the MFA the weak 
FM produces a divergence of the z-component of the sus- 
ceptibility only in a very narrow temperature region close 
to the Nccl point. Since the MFA does not account for 
near-neighbour correlations between the spins, away from 
the immediate vicinity of Tm the weak FM is ignored and 
behaves like a T-independent transverse susceptibil- 
ity (that is, transverse to the ordered AF moment). In 
contrast to this, the z-component of the susceptibility 
calculated within the RPA has a strong temperature de- 
pendence and shows that the effects of the quantum fluc- 
tuations are important in a wide region below the Neel 
temperature. The other component which shows some 
differences between the MFA and RPA below the transi- 
tion is ^ Slid for this component it is seen that since the 
MFA does not include the reduction of the staggered mo- 
ment (which is in this direction) due to quantum fluctu- 
ations at low temperatures, linear spin- wave theory, and 
not the MFA, agrees with the RPA for low temperatures. 

Further, above the Neel temperature the differences 
can be understood as follows. In a MFA (that is T^^^^ « 
J), both components of the susceptibility x^ and x^ are 
rapidly changing functions in the immediate vicinity of 
Tjv {5T/Tm ~ 0.005), and then have the same behaviour 
as the term further above the transition. This MFA 
behaviour is in no way similar to that found in the RPA. 
That is, our results are an example of the pronounced 
effects of short-range correlations and quantum fluctua- 
tions. The RPA scheme gives a much lower value of the 
Neel temperature than MFA does (Tn ~ 0.3 J), but in a 
broad T region above the Neel point strong short-range 
correlations exist, and the RPA includes the manner in 
which these fluctuations strongly modify the suscepti- 
bility. Similar reasoning explains the differences in x^ 
between the MFA and RPA. 

For completeness, in Figs. llUlllll2l we present all com- 
ponents of the susceptibility together, within both the 
MFA and the RPA, contrasting different values of the 
physical parameters describing the DM interaction. To 
be specific, in Fig. ^| we show the situation when ei > 62 
with the ratio (ei/£2)^ ~ 4.2 at zero temperature. As a 
result, we obtain that for T = < X^iX^ with the 
same ratio between the x, z and y components of sus- 
ceptibility x'^'Vx^ ~ 4.2 (see flV^ . By increasing the 
magnitude of the DM parameter d, due to the strong 
dependence of the mode £2 of the d (see Fig. 0J|, we ob- 
tain the situation corresponding to £2 — £\- Then, as is 
seen in Fig. ^] both the MFA and RPA schemes result 
in equal values of the all components of susceptibility at 
low T. A further increasing old leads to the situation 
£2 > £1, opposite to the one presented in Fig. ^1 In 
Fig. El we show the susceptibility in the case of the ra- 
tio (£2/£i)^ « 2.0, and at T = one finds x^ > X'.x'' 
and x^ Ix^'^ ~ 2.0 for T = 0. (We note that for other 
sets of d, AF, the behaviour of the components of x is 
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FIG. 10: (Color online) AH 3 components of the susceptibility 
within (a) the MFA, and (b) within the RPA, for d/J = 0.02 
and Ar/J = 0.42 x 10"^ 



determined almost entirely by the ratio of the spin- wave 
gaps, £i/£2 =. These results agree with our analytical 
predictions (see section IV).) 

Then, comparing the z-component of the susceptibil- 
ity within the RPA (Figs. [Tni;b)-|T2;b)) we also find that 
increasing of the anisotropy parameter d leads to the 
broadening T regions where (i) the effects of the quantum 
fluctuations are important T < Tn and (ii) the strong 
short-range correlations exists T > Tn. 



VI. SUMMARY AND 
CONCLUSIONS/DISCUSSION 

To summarize, we have presented a theoretical in- 
vestigation of a single Cu02 plane of the undoped 
La2Cu04 crystal in the low-T orthorhombic phase. The 
Cu spins in the plane were modelled by the 2D spin- 
1/2 Heisenberg AF with spin-orbit coupling, the lat- 
ter represented the antisymmetric and symmetric DM 
anisotropies. We have adopted the Green's function 
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FIG. 11: (Color online) All 3 components of the susceptibility 
within (a) the MFA, and (b) within the RPA, for d/J = 0.041 
and Ar/J = 0.42 x 10"^ 
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FIG. 12: (Color online) All 3 components of the susceptibility 
within (a) the MFA, and (b) within the RPA, for d/J = 0.058 
and AF/J = 0.42 x 10"^ 



method within Tyablikov's RPA decoupling scheme to 
calculate the magnetic susceptibility of such a model. 
In order to allow us to accurately model the longitu- 
dinal susceptibility within such a level of decoupling of 
high-order Green's functions, we have extended Liu's 
methodic for the isotropic Heisenberg model to one that 
includes a weak canted FM moment in the plane. 

We can emphasize several important conclusions from 
our results. We have found that the anisotropy in- 
troduced into the problem by the symmetric and anti- 
symmetric DM interactions leads to important changes 
in the behaviour of the magnetic susceptibility near the 
transition point. By comparing the MFA and RPA re- 
sults we conclude that the effects of quantum fluctua- 
tions and the short-range correlations are very strong in 
the broad temperature region of near the Neel tempera- 
ture. Further, we find that since the RPA and SW re- 
sults are quite different near the Neel temperature, the 
effects of spin-wave interactions, which are included in 
an approximate way in the RPA but not the SW the- 
ories, are very important in this system. This neces- 



sarily leads to the question, would more advanced de- 
coupling schemes, namely improvements on Tyablikov's 
decoupling {e.g., see our Eq. l|5^ ). or, possibly, the inclu- 
sion of nonlinear effects in the SW theory, lead to quali- 
tatively different results? 

Secondly, we have obtained that the weak fcrromag- 
netism in the z-direction (caused by the DM interaction) 
leads to the essential difference between the temperature 
behaviours of the transverse s-nd components of 
the susceptibility (recall that the AF moments lie in the 
y — z plane and are nearly aligned along the y axis). We 
established the correlation between the ratio of the in- 
and out-of-plane spin-wave modes of the excitation spec- 
trum in the long wavelength limit (fc = 0), which is fixed 
by the ratio between the d and Pi — Pa DM parameters, 
and the behaviour of x^''^ ''■s. x^ ™- the zero temperature 
limit. This conclusion is independent on the analytical 
method which we used to calculate the susceptibilities, 
since all methods agree in the low-T regime, and could 
allow one to make predictions concerning the gaps in the 
excitation spectrum based on the data for the suscepti- 
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bility. 

Now we comment on the comparison of our results to 
the experimentally observed anisotropies^ that motivated 
this work. We can state that, in addition to the known 
resultsiSiiii^ that DM interaction induces the weak fer- 
romagnetism in the LTP phase and the spin- wave gaps, 
this interaction is at least in part responsible for the un- 
usual anisotropy in the magnetic susceptibilityiiS, We can 
mention the most significant features observed in the ex- 
periment that are in qualitative agreement with the pre- 
sented in paper theoretical results: (i) the absence of any 
special behaviour (anomaly) in the transverse component 
across the Neel temperature; (ii) the additional in- 
crease of the x** component in the ordered state and its 
smooth decrease in a broad temperature region in the 
paramagnetic state; (iii) a significant temperature de- 
pendence of the component in the broad temperature 
region below and above the transition point. 

Now we briefly discuss the experimental data which 
cannot be explained within the framework of the pro- 
posed here theory. Firstly, we have found that the ob- 
served ratio between the x and y components < 
(in the T = limit) takes place only if the spin-wave 
gap with out-of-plane mode is less than the in-plane one 
Eo < £i. However, older neutron-scattering experiments^ 
find the opposite ratio: the magnitude for the out-of- 
plane mode is 5 meV, for the out-of-plane mode 2.3 meV. 
Recent Raman work confirms one of these valuesi^ So, 
other interactions which affect these gaps must be im- 
portant for an accurate explanation of the susceptibil- 
ity data. Secondly, our results cannot explain a T- 



independent shift between x^ , X^ ^-i^d x^ observed in ex- 
periments - an explanation of this physics is provided in 
the experimental paper, namely that one must include a 
van Vleck contribution which shifts, in a T-independent 
manner, these components of the susceptibility, but we 
defer our inclusion of this physics until the second paper 
in this series of theoretical studies. 

For further improvements of our theoretical modelling 
of the La2Cu04 compound, it seems to be important to 
investigate a 3D model on a body-centered lattice with 
the weak AF interlayer coupling. It is also possible to ex- 
tend the 2D model by considering the ring exchange, and 
the interaction between the next nearest neighbour sites, 
and we expect that some of these additional physics can 
be responsible for the correct ratio between the spin- wave 
gaps with respect to the ratio between x^ and . In ad- 
dition, the anisotropic Van Vleck contribution (orbital 
susceptibility) and gyromagnetic (Lande) factor need to 
be taken into account. We will present a detailed compar- 
ison to these experiments when these other interactions 
are included in future publications. 
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